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Abstract 

A new global two-fluid electromagnetic turbulence code, CENTORI, has been developed for the 
O-l purpose of studying magnetically-confined fusion plasmas on energy confinement timescales. 

This code is used to evolve the combined system of electron and ion fluid equations and Maxwell 
equations in toroidal configurations with axisymmetric equilibria. Uniquely, the equilibrium is 
co-evolved with the turbulence, and is thus modified by it. CENTORI is applicable to tokamaks of 
arbitrary aspect ratio and high plasma beta. A predictor-corrector, semi-implicit finite difference 
scheme is used to compute the time evolution of fluid quantities and fields. Vector operations 
and the evaluation of flux surface averages are speeded up by choosing the Jacobian of the trans- 
formation from laboratory to plasma coordinates to be a function of the equilibrium poloidal 
magnetic flux. A subroutine, GRASS, is used to co-evolve the plasma equilibrium by comput- 
ing the steady-state solutions of a diffusion equation with a pseudo-time derivative. The code is 
written in Fortran 95 and is efficiently parallelized using Message Passing Interface (MPI). Il- 
lustrative examples of output from simulations of a tearing mode in a large aspect ratio tokamak 
plasma and of turbulence in an elongated conventional aspect ratio tokamak plasma are provided. 

Keywords: Two-fluid and multi-fluid plasmas, Drift waves, Tokamaks, spherical tokamaks, 
Plasma turbulence, Magnetohydrodynamic and fluid equation 

PACS: 52.30.Ex, 52.35.Kt, 52.35.Ra, 52.55.Fa, 52.65.Kj 

1. Introduction 

Plasma confinement in tokamak experiments is determined partly by binary Coulomb colli- 
sions between charged particles, but mainly by turbulence and instabilities, which occur on scales 
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ranging from particle Larmor radii to the system size. Understanding the nature of this turbu- 
lence is a key goal of thermonuclear fusion research, since the confinement time is one of the 
parameters that must be optimised in order to create burning plasma conditions. In order to sim- 
ulate turbulence in tokamak plasmas it is necessary to either average the Vlasov equations of the 
particle species over gyro-angle (the gyrokinetic approach) or take full velocity-space moments 
of these equations (the fluid approach). The lower dimensionality of fluid models makes it possi- 
ble to simulate larger systems over longer timescales, and for this reason fluid codes continue to 
play an important role in tokamak plasma modelling. Some of these codes are based on electro- 
static models 01 [2] or employ flux tube geometry [3], while others are designed specifically for 
the purpose of simulating edge plasma phenomena, such as edge localised modes (ELMs) |4, 5]. 
A global magnetohydrodynamic (MHD) code NIMROD [6] has also been applied to the mod- 
elling of ELMsj7[], in addition to a range of other MHD instabilities in several different toroidal 
configurations |8y. In order to model turbulent transport on confinement and resistive diffusion 
timescales in an electromagnetic global code, it is necessary to include two-fluid effects, and it is 
also desirable to co-evolve the equilibrium. 

In this paper we describe CENTORI (Culham Emulator of Numerical TORI), a new toroidal 
two-fluid, electromagnetic turbulence simulation code that meets these requirements. It can be 
used to describe the co-evolution of turbulence, MHD instabilities and equilibrium in tokamak 
plasmas with arbitrary aspect ratio and high plasma beta (ratio of plasma pressure to magnetic 
field energy density). It is designed for the specific purpose of simulating global two-fluid elec- 
tromagnetic tokamak plasma turbulence on confinement timescales, in realistic geom etries and 
in conditions such as those found in the present-day machines MAST |9] and JET II X Oil - and in the 
forthcoming international fusion experiment ITER 111 ill . Turbulent modes in tokamak plasmas 
are typically drift waves, which are predominantly electrostatic waves driven by temperature or 
density gradients. An important example is the ion temperature gradient mode, which has wave- 
lengths perpendicular to the magnetic field of the order of the ion Larmor radius p, 1 12]. Many 
tokamak turbulence codes, such as the electrostatic fluid codes mentioned above and also gyro- 
kinetic codes such as Kinezero lfl3ll . are designed specifically for the modelling of drift waves 
in a fixed, prescribed plasma equilibrium. CENTORI, on the other hand, is designed to study the 
interaction between drift waves and MHD instabilities, which generally occur at longer wave- 
lengths, ranging up to the system size, in a co-evolving equilibrium. However fluid codes such 
as CENTORI cannot be used to model explicitly instabilities that occur on the smallest tokamak- 
relevant spatial scales, in particular length scales below the ion Larmor radius. Phenomena on 
the scale of the electron skin depth S e are specifically excluded from the model used in CENTORI, 
since electron inertia is neglected (in any event 6 e < pi unless the plasma beta is less than the 
electron to ion mass ratio, which is not normally the case in the core region of tokamak plasmas). 
The drift waves described by gyro-kinetic theory have frequencies of the order of p*Q where 
p* is particle Larmor radius normalised to the equilibrium gradient scale length and Q. is the 
corresponding cyclotron frequency lfl4ll . Two-fluid theory, on the other hand, can accommodate 
MHD modes such as global Alfven eigenmodes fl5ll . which, in low beta plasmas, have frequen- 
cies higher than those of ion drift waves. CENTORI can be used to study processes occurring on 
timescales ranging from the reciprocal Alfven frequency to the energy confinement time. 

The physics model implemented in CENTORI is very similar to that used in CUTIE, a global 
two-fluid electromagnetic turbulence code which was based on periodic cylinder geometry and 
was restricted to large aspect ratio plasmas with circular poloidal cross-section 11611 . Despite 
these restrictions, CUTIE has been used for a number of successful applications. For example, 
it was recently shown to reproduce experimentally-observed transitions to a high confinement 
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mode of plasma operation via the control of particle fuelling in the COMPASS-D tokamak lfl7ll . 

This paper is organised as follows. In Section|2]we describe the relationship between labora- 
tory coordinates and plasma coordinates, in which the fluid and Maxwell equations are evolved 
in CENTORI. The form in which these equations are solved is discussed in Sections [3H5J while 
initial and boundary conditions are discussed in Section [6] Sections [7] and [8] are concerned re- 
spectively with the distinction made in the code between mean and fluctuating quantities, and 
global quantities evolved by it, such as plasma beta. In Section [9] we describe GRASS, a sub- 
routine of CENTORI which co-evolves the plasma equilibrium using a novel pseudo-transient 
method. Operational and technical aspects of the CENTORI package and the code structure are 
discussed in Section [10] while in Section QT| we present some illustrative examples of output 
from a simulation of a large aspect ratio tokamak plasma. 

2. Coordinate system 

Before describing the physical quantities and their evolution equations, it is useful to provide 
a full description of the coordinate systems used in CENTORI. 

2.1. Laboratory coordinates 

CENTORI is used to model a toroidal plasma held in place by magnetic fields produced by 
external coils and by the plasma itself. A natural coordinate system to use for the laboratory 
frame is the right-handed cylindrical system (R, Z, f), where R is major radius (distance from the 
machine's vertical axis of symmetry), Z is vertical distance (parallel to the symmetry axis), and 
£ is toroidal angle (azimuthal angle around the symmetry axis). We note that 

V<r = -V0 = ie f , (1) 

where <p is azimuthal angle in the right-handed cylindrical system (R, <f>, Z) and is the unit 
vector in the ( direction. 

2.2. Plasma coordinates 

The total magnetic field in the system comprises the vacuum field, produced solely by cur- 
rents flowing in conductors surrounding the plasma, plus the field generated by the currents in the 
plasma itself. We use the total equilibrium magnetic field to define the plasma coordinate system. 
The equilibrium poloidal flux function if/(R, Z) defines the equilibrium poloidal magnetic field. 
The quantity t//(R, Z) can evolve in a CENTORI simulation, but only on a much longer timescale 
than the turbulence. It is the magnetic flux per unit toroidal angle passing through the horizontal 
circle of radius R centred at (R = 0,Z); it is independent of When plotted in the poloidal 
(R, Z) plane the lines of constant if/ in the vicinity of the plasma form nested, closed contours 
(flux surfaces). The minimum value of i// within these closed surfaces lies near the centre of the 
plasma, and defines the location of the magnetic axis, along the circle (Rq, Zo, £). 

In a real machine the edge location of the plasma is determined by either a physical limiter or 
the design of the magnetic geometry. Because only the gradients of if/ have physical meaning we 
may, for convenience, adjust if/ so that the known location of the edge of the plasma is defined 
to lie on the if/ = contour. Figure Q] shows a typical set of if/ contours in the poloidal plane. 
This plot is effectively the starting point for the calculations performed using CENTORI. The flux 
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Figure 1 : Typical plot of if/ contours over the (R, Z) grid, showing laboratory coordinates (R, Z, Q and plasma coordinates 
(if/, 0, f ) employed in CENTORI. This plot was obtained using the GRASS equilibrium solver (Section f9.lt . 



contours if/(R, Z) are determined a priori either by an external program or the equilibrium solver 
in the code, which is described in Section |9j they are co-evolved in time with the turbulence. 

Our aim is to evolve a set of plasma quantities which are stored in arrays at a convenient set of 
computational grid points in a right-handed but in general non-orthogonal dimensionless plasma 
coordinate system (p, 9, £). Here p is a radial coordinate, with Vp directed from the magnetic 
axis to the plasma edge, and 9 denotes an angle in the (R, Z) plane. 

2.3. Radial coordinate 

The radial coordinate p is a normalised measure of iff, the normalising factor being the abso- 
lute value of the poloidal flux at the magnetic axis, iffg. Thus, from the magnetic axis at (Ro,Zo) 
to the edge of the plasma we have -iffo < iff < and < p < 1 with p defined in terms of ift by 

P = 1 + iA/«Ao ■ (2) 

The radial grid points are equally spaced in p. In the cylindrical limit p varies approximately as 
r 2 where r is distance from the magnetic axis. The p contours are thus relatively far apart near 
the magnetic axis, as shown in Fig.Q] Because the magnetic axis is a coordinate singularity (all 
9 points at p = and a given ( coincide), we have chosen to locate the innermost p grid points 
on a contour that is slightly displaced from the axis itself. 

The gradient Vp in the laboratory frame is determined from the ifr(R, Z) grid by fitting two- 



dimensional Chebyshev polynomials [18] to the known iff values at the grid points, and taking 



their derivatives in the R and Z directions. It follows from Eq. (f2j) that 



Vp= — V<A = — ^e R + — -|e z , (3) 
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where eR and ez denote unit vectors in the R and Z directions. 

2.4. Relationship between equilibrium magnetic field and plasma coordinates 
The equilibrium poloidal magnetic field is given by 

B p = xV^ = ^ (Vf xVp). (4) 

Thus = iffo\Vp\/R. The toroidal equilibrium magnetic field is given by 

B, = FVf, (5) 

where the scalar quantity F is taken to be a flux function, i.e. it depends only on the radial 
coordinate p. This is generally a good approximation under typical tokamak conditions lfl9ll . 
Thus the total equilibrium magnetic field is 

B eq = ^ (Vf xVp) + FV£. (6) 

We define a vector potential A in the usual way as a vector field whose curl is equal to the 
magnetic field. We can write the equilibrium vector potential A eq in covariant form as follows: 

A eq = A eqp Vp + A eq e V6» + A eq ( Vf . (7) 

For convenience we choose a gauge such that the radial component of A eq vanishes, i.e. 

Aeq P = 0. (8) 

In terms of the remaining components of A eq , the equilibrium magnetic field becomes 

B eq = VxA eq = VAeq e xV0 + VA e q f xV(. (9) 

Matching the poloidal components of Eqs. (0 and (0 we find that we can set 

Aeq^ = -i/r. (10) 

Matching the toroidal components of Eqs. © and (0 we obtain 

dAeqe 

FV£ = VAeqsX V0 = — — VpxV0, (11) 
n op 

and the scalar product of this with V£ yields 

F dA e qg dA e qg 
FV(-V{= — = —^-Vt-(VpxV6) = —^-J, (12) 

where J e V( ■ (Vp x V0) is the Jacobian relating laboratory and plasma coordinates (see 
following subsection). The covariant poloidal component of A eq is thus given by 

Aeq fl = J j^dp. (13) 

Eqs. <j8j, (TTOl i and (fT3l define the equilibrium vector potential A eq in covariant form; the equilib- 
rium magnetic field B eq may be calculated by taking its curl. 

The set of space variables (p, 9, £) constitutes a quasi-orthogonal coordinate system in which 
Vp ■ V£ = V6> • V^" = 0, but in general Vp • V# + 0. Taking scalar products of B eq with the 
coordinate gradients we obtain 

B eq -Vp = 0, B eq -V# = (Ao(V^xVp)-V0 = iAoJ", B eq ■ V£ = F/R 2 . (14) 

These three equations give the contravariant components of B eq directly, thereby eliminating the 
need to perform a curl operation (see Section l23i l. 
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2.5. Poloidal coordinate 

The poloidal angle 9 varies from to 2n in the (R, Z) plane. By convention, points at 9 - 
lie along the line defined by (R > Rq,Z = Zq), and 9 increases in the anticlockwise direction as 
shown in Fig.Q] Denoting by I the arc length in the poloidal plane along a given p contour, we 
can write 

89 IVpl 89 

B ai -W = B P — = =^J. (15) 

To determine the distribution of 9 grid points along the contour in the (R,Z) plane we introduce 
a parameter r and solve the following pair of Hamiltonian equations 1 1911 : 

dR dp dZ dp 
dr dZ ' dr dR' 

with (R,Z) being stored at intermediate points as the solution proceeds. The gradients in p are 
calculated using Chebyshev polynomials, as described above, and a convergence loop ensures 
that the contour is followed with sufficient accuracy. The arc length I is given in terms of r by 



We choose J' to be a flux function, i.e. J' = J'(p). This enables 9 points on a given flux contour 
to be determined by integrating the expression 

d9 = RJ-^-=RJdr, (18) 
|vp| 

where J' is obtained by imposing a 2n periodicity on 9: 

Jip) = < 19) 
§RdT 

It is straightforward to interpolate the stored (R, Z) values to determine the locations of equally- 
spaced 9 points along the p contour. Figure [2] shows an example of a (p, 9) grid. 

The process described above can be used to map out the locations R(p, 9), Z(p, 9) along the p 
contours. The partial derivatives dR/dp, dR/d9, dZ/dp and dZ/d9 are found by fitting Chebyshev 
polynomials to R and Z along the p direction and Fourier series in the 9 direction. These provide 
the contravariant basis vectors of the plasma coordinate system: 

b p = T (ye xV^ = ^e R +^ e z , (20) 
dp dp 

b fl = J"*(V^xVp) = ^eR + ^ez, (21) 
d9 d9 

b ? =J"(VpxV0)=Re o (22) 

where J'* = 1 1 J' and follows directly from Eq. ([]]); V£ is the covariant £ basis vector and 
hence is reciprocal to b^ . We may then calculate J'* (and therefore JJ) using 



J* = be ■ (b f x b p ), 
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(23) 




Figure 2: Typical set of (p, 6) grid points, superimposed on the original ip(R,Z) grid. In this case there are 129 radial grid 
points and 65 poloidal grid points. 



and the covariant basis vectors are given by 

bP = Vp=J(b e xb ( ), b e = V6=J(b ( xb p ), tfsVf=i^. (24) 

It is straightforward to evaluate these vector products since the covariant and contravariant basis 
vectors are all stored with components in the laboratory frame (although they are evaluated at 
specified points in (p, 6) space, their components relative to the basis (eR, ez, e^) are known). 

We adopt this particular algorithm to obtain the gradients and the Jacobian in order to max- 
imise accuracy and smoothness in the results through the use of the Chebyshev/Fourier fitting 
method, and it also guarantees that the covariant and contravariant basis vectors are reciprocal. 

2.6. Vector operations in plasma coordinate system 

This section provides expressions for scalar and vector products together with differential op- 
erators in the plasma coordinate system. In what follows A and B are arbitrary vector functions, 
while / is an arbitrary scalar function. The vector A has covariant representation 

A = A p b p +A fl b fl +A f b f , (25) 

where A, ■ = A • b, . The corresponding contravariant representation is 

A = A p b p +A e b +A i; b ( , (26) 

where A' = A • b'. The scalar product of A and B is then A ■ B = A,B' = A'B,, where a repeated 
index implies summation, while the vector product is given by 

A x B = J^A e B ( -A ( Bg)b p + (A ( B p -A p B ( )be + (A p B e -A e B p )b^ 

= J* {(A 9 B ( - A^B e )b p + (A ( B p - A p B ( )b + (A p B e - A B p )b ( \ . (27) 



The gradient operator, which is defined in the usual way, produces a covariant vector. The diver- 
gence of A is evaluated using its contravariant components: 



while the curl is obtained using its covariant components, and the result is a contravariant vector: 



u dA, dA g \ ldA p OA A IdAg dA p , . 



The choice of J' as a flux function considerably simplifies and speeds up many calculations. 

2. 7. Physical coordinates 

It is convenient to perform the vector operations discussed above using the covariant and 
contravariant representations. However, these do not have directions, dimensions or units that 
are intuitive as far as the physics is concerned. We therefore define a third set of components for 
the vector quantities in CENTDRI, which we refer to as their physical representation: "normal", 
denoting the direction normal to the flux surface; "tangential", parallel to the flux surface in the 
(R,Z) plane; and "toroidal", around the machine axis. The physical components are orthogonal: 

Vp _ A p 

A„,=A.— = — , (30) 
(V<TxVp) (V^xVp) _ JR 

Atangential -A- ^^-A- - Ag ^ , (31) 

A _» I_f ™ 

Atoroidal - A ■ ^ - ^ . 

2.8. Flux surface-averaged quantities 

It is necessary to compute flux surface-averaged quantities in CENTDRI since these affect the 
evolving equilibrium. The flux surface average of a scalar quantity f(p, 9, g ) is given by 



JOT f do da j 
'' ~ ff dodci:! 



a >> = "rr ... . = ^ JJ f, '"'' c (33) 

where / is evaluated at fixed p and we have used the fact that J is defined to be a flux function. 



3. Physics quantities 

The primary quantities evolved by CENTORI are as follows: Vi, ion velocity; A, vector po- 
tential; rii, ion number density (= n e , electron number density, via quasi-neutrality); T,, ion 
temperature; and T e , electron temperature. In addition, a number of auxiliary quantities can be 
advanced in time once the primary quantities have been updated. These are: v e , electron veloc- 
ity; <D, electric potential; B, magnetic field; J, current density; E, electric field; p, = n e Tj, ion 
pressure; and p e = n e T e , electron pressure. 
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These variables are normalised as follows: 

v^, v: = ^, A* = A, b-=» (34) 
va v a do "o 

n e = =, Ti = — , T e = —, Pi = n e T i = — , p e = n e T e = — , (35) 

n e / 10 *e0 PiO PeO 

where va = Bo/ -\j4np m Bo/ ^<Xtt ra,- ~ e is a typical Alfven speed, p m = m, n e is the ion mass 
density, Z?o is the vacuum toroidal field at the magnetic axis, ~ e is the volume-averaged electron 
number density, r,o is the initial ion temperature at the magnetic axis, T e o is the initial electron 
temperature at the magnetic axis, p/o - 7T e Tg) is a nominal ion pressure, p e o = 7T e T e o is a nominal 
electron pressure, and m, is the ion mass. The quantities Bo, r,o and T e o are given nominal values 
by the user; ~ e is calculated as the plasma evolution progresses, so va, pio and p e o vary with time. 
It should be noted that the actual density and temperature values on axis are not constrained to 
be their initial arbitrary values, but vary as the profiles evolve. The normalised quantities listed 
above are all dimensionless except for A* which has the dimensions of length. 



4. Two-fluid equations 

4.1. Momentum equations 

The ion momentum balance equation can be written in the form 

Pm (j£ + W x ViJ = -V Pi - ^Vvj 2 + en e E + ^(v, x B) - en.rH - p mX vF x W) + S v . (36) 

where W = V x Vi is vorticity, rj is resistivity (assumed to be a scalar function of space and time), 
Xv is velocity diffusivity (see Section 15.1.21 ). S v is external force density (see Section 15.1.21 ). e 
is proton charge and c is the speed of light (we use Gaussian cgs units throughout this paper, 
although output from the code is in SI units, to facilitate comparison with experimental results). 
In the electron momentum balance equation we neglect inertial terms, momentum sources and 
viscosity: 

= -Vp e - en e E ^(v e x B) + en e r]J. (37) 

c 

This is equivalent to Ohm's law in the limit of vanishing electron mass. 

4.2. Energy equations 

The transfer of energy is described by the two equations 

\ Ue [it + Vi ' v ) Ti + PiV ' Vi = ~ v ' q '' + 5 '* (38) 

\n e [j t + v e • vj T e + Pe V ■ v e = -V • q e + S e , (39) 

where q,; e are the ion and electron heat fluxes (see Section 15.41 ) and S ,> are additional ion and 
electron heating sources. 
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4.3. Mass continuity equation 

The mass continuity equation used in CENTORI is 



dp m 
dt 



+ V ■ (p m Vi) = S n - nij n e v A V -T* w + S„ - v m (p m - <p m », (40) 



where S„ is the particle source rate (see Section l5.5.U . V • T^, is a term representing the effect of 
the Ware pinch [20] (see Section l54l i. 6„ is a diffusion term given by 

with^„ e and^RR respectively the particle and Rechester-Rosenbluth diffusivities (see Sections [5Tl 
and !5.5.2l i. Finally in Eq. d40l l. v/y is the parallel ion thermal relaxation rate (see Section l5~4l i. 
It is not necessary to solve an electron continuity equation since the plasma is assumed to be 
quasi-neutral and the current is assumed to be divergence-free. 

4.4. Maxwell's equations 

The vanishing of the divergence of B is guaranteed in CENTORI through the use of the poten- 
tial representation B = V x A and the induction equation is solved for A rather than B: 

1 <9A 

-— =-E-VO. (42) 
c dt 

Current densities J are computed using the pre-Maxwell form of Ampere's law: 

J-fvxB. (43) 

47T 



5. Normalised physics equations and their solution 

CENTORI is used to evolve the normalised quantities defined by Eqs. d34l l and (l35l l rather 
than the absolute values of velocity, magnetic field, and so on. In this section we explain how 
the physics equations themselves are normalised. Unless otherwise stated, all of the normalised 
equations have the dimensions of reciprocal length. The equations are solved by using finite dif- 
ferences to approximate all of the derivatives; the solution method is thus entirely non-spectral. 
A key advantage of this approach is that parallelisation of the code is then relatively straightfor- 
ward, and yields good scalability results (see Section ITOb. On the other hand the finite-element 
method, used, for example, in NIMROD |@], is particularly well-suited to modelling the edge re- 
gions of plasmas with strongly-shaped poloidal cross-sections. 

5.1. Normalised ion momentum equation 

All three physical components of \* are evolved, with subscript "1" labelling the normal 
direction, subscript "2" the tangential direction, and subscript "3" the toroidal direction. We 
define a normalised vorticity W*: 

W 

W* = — = (V x v.* ). (44) 
va 
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It should be noted that W* has the dimensions of reciprocal length. Using also the normalizations 
introduced previously, dividing by v\m{fQn*, defining the following quantities 



B i0 = 



4tt pjo 
Bi 



— , PeO = 



v\ m -, n. 



Bi 



TeQ 

v\m,i 



n _ Xv c « _ S v _ 4ttS 

v A v a nii n e Oq 



and introducing an additional term related to the Rechester-Rosenbluth diffusivity Drr 112 111 [see 
Eq. (17 1H . we find that the ion momentum equation can be written in the form 



±ch± 

va dt 



va 



1. 



S. 



«* i n* 



va q \ (K) 



v<«:> 



D v (V x W) - 



va 



1 <9A* 



(45) 



where J* = 47rJ/(cZ?o), if = c 2 t]/(47tva), and co C i = eBo/QtiiC). To reduce problems arising 
from short wavelength modes in the radial direction, the momentum equation is supplemented 
by artificial damping terms: 



v A dt I i, normal j *" " i, 



i. normal' 



(46) 



where 6 V = 0.5Vj-|[/va, Vj-y being the parallel ion thermal relaxation rate [Eq. (l75lll. A similar 
damping term is applied in the tangential direction. The dimensionless damping rate used in the 
code is SI = v^Afc^,. 



5.1.1. Evolution of normalised momentum equation 

In the current version of CENTORI we neglect the (1/va)0A* /dt and J* terms on the right 
hand side of Eq. d45l ). In tokamak plasmas there is generally a large separation between drift and 
Alfven timescales, with the consequence that turbulent fluctuations are predominantly electro- 
static in nature, and the inductive part of the electric field term plays only a minor role in the ion 
momentum equation. In Section [Tl.2l we will illustrate this point using results from a CENTORI 
simulation. The J* term in Eq. d45b is small by virtue of the fact that tokamak plasmas tend to 
have very high Lundquist numbers, i.e. are close to being perfectly conducting. 

In Eq. d45b it is not straightforward to convert 6 = -£\,(V X W*) into finite differences, due 
to the non-orthogonal nature of the coordinate system. We approximate it by the expression 



6*R JD r 



/ 5 2 v* {B p ) 2 d 2 y{ 
~dp T + Bj ~W 



(47) 



We are assuming here that the contribution of viscosity to the ion momentum equation can be 
well-approximated by a term proportional to V 2 v*. We are thus neglecting the V(V ■ v*) term 
in V x W* (although the flows described by CENTORI are in general compressible); it is not 
necessary to include this term in order to model the neoclassical and turbulent damping of flows 
ll22ll . The factor containing B p is present to take account of the spacing of adjacent points in the 
6 direction being proportional to the local poloidal field [see Eq. dl5H . Adopting the convention 
that subscripts /', k and I label array elements in the p, 9 and ( directions respectively, while 
superscripts label time, we approximate 6 by the central difference expression 



^0 Jj,k Dv j.k,l 

(ApJ 2 



ft 



N+l *N+\ 
;'+l,jy + V i j-\,k,l ' 
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■2y* N+l \ 
zv ij,kj) 



J<k(u l>> a. (Bp) 2 / N N _ 2 *n + i\ (4g) 

(AO) 2 B 2 V'jMU + y ij,k-U ^ijXi)- ^ 

We also introduce purely numerical diffusion terms, with coefficients e p , eg and e^, which are de- 
signed to remove variations in v? of similar length scale to the separation of adjacent grid points. 
It is evident that the resultant finite-difference equations are consistent with the governing partial 
differential equations as the mesh sizes tend to infinity. These effectively suppress spurious os- 
cillations at wave numbers corresponding to inverse mesh size. Unlike the turbulent diffusivities, 
they are non-zero even when the turbulent fluctuation amplitudes go to zero for any fixed mesh 
size. Thus the finite difference approximation to the momentum equation is of the form 

Here, superscripts L ("latest") indicate the most up-to-date (most time-advanced) values avail- 
able. This is to avoid the use of "new" values at adjacent 9 and f points (i.e. at k± 1, 1 ± 1); these 
would appear as undesirable off-diagonal terms in the tridiagonal matrix equation. Typically, 
\* L = v. from the previous iteration. The numerical diffusion coefficients e p , eg and e f have 
the following forms: 

y/p _ 1 

€p ~ N 2 ' €8 ~2tt 2 N 2, e( ~87T 2 N 2 ' 

where N^, Ng and N( are the numbers of grid intervals in the respective directions. 

Dropping the k and / subscripts on velocity components, the finite difference approximation 
to the momentum equation can be written in the block tridiagonal matrix equation form 

A^+B^+C^-Rj, (50) 

where A ., B ., C . are 3 x 3 matrices and Rj is a vector that depends on the latest (L) values of 
velocity components as well as those of the previous timestep. Equation d50l > is solved for the 
normalised ion velocity vf at the new timestep using a standard predictor-corrector scheme, with 
L . converging to vj* . . We then subtract the flux surface-averaged normal component of Vj , so 
that only a fluctuating part remains. 

5.1.2. Momentum sources and transport 

Presently only a toroidal momentum source is included in CENTORI; the profile is given by 

r-* / \ _ f 4^ Paux,i(p) + Paux,e(p) 

" vjoAP) = J mom ~^ T^TT i P-l) 

Bq v th j(0) 

where P m u.i/e is the external heating power per unit volume provided to the ions/electrons, 
VtA,i(0) = (2r,(0)//«,) 1 '' 2 is the ion thermal speed at the magnetic axis and f mom is a user-defined 
multiplier which matches the total momentum provided to the plasma with experiment. 

Turbulent diffusivity terms are included in the full definition of the normalised velocity dif- 
fusivity D v [see Eq. d45H : 

D v (p,0,O = *22= ( 1 + q(R} 2 ,&/yy J*' + W* 2 l) + ^±ffi££ ( 52 ) 
v A \ \m e l 1 v A 



- 2 ~ 2 

where ;^,, iUser is user-specified, J* is a normalised measure of entropy and W*" is a normalised 
enstrophy (J* and W* being the fluctuating parts of the normalised current density and vorticity 
respectively). The parameter fjj is a user-defined multiplier between and 1. The final term in 
D v is a Gyro-Bohm diffusivity: 

( ^ t V,Kip2 i f ™ 

^v.classicaltP) - Jyc . (JJ) 

a 

with < f xc < 1 a user-defined multiplier and p, the ion gyroradius. 

5.2. Evolution of normalised electron velocity 

The normalised electron velocity v* is determined directly from \* and J* by noting that the 
net current density J is given by 

J = e n e (Vi - v e ). (54) 

Hence 

v; = v?- c l\ r. (55) 

47T e n e n* v& 

5.3. Evolution of electromagnetic quantities 

5.3.1. Normalised Ampere 's Law 

It is evident from the definitions of B* and J* that Ampere's law has the normalised form 

J* = VxB*. (56) 
Note that B* is dimensionless while J* has the dimensions of reciprocal length. 

5.3.2. Normalised Faraday law and Ohm 's law 
Ohm's law [Eq. ( 1371 )1 can be written in the form 

V e XB V Pe 

E = + 77J. (57) 

c en e 

We divide the electric field into ideal and resistive parts by writing 

_ V e XB V Pe 
•iideal , 

c en e 
E res = Tli, 

and write Faraday's law in the form 

- -T- = -Eideal - V<D - E res = - - ^ £i - V0> - 7/J. 

c ot \ c en e 



We can thus write 



c 



dt \ c e e n e n* I An 
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and hence, multiplying by c/(v,a Bq), we obtain 



1 dA 

va dt 



v* x B* + 



v A Boen~ e n* 
We thus obtain the normalised Faraday's law 

1 dA* „„ cT e o 



v A dt 



-E* 



v a Bo e 



vaBq 



Anv J 



■nT. 



(58) 



where the normalised electric potential <f>* is defined as e(&/T e o and the normalised electric field 
E* is defined by a normalised Ohm's law 



E* = 



va Bo 



-E 



-v* x B* 



C TeQ Vp* e 



v A B e n* 



+ n*r, 



(59) 



the term in parentheses being the ideal part and the remainder the resistive part. 



5.3.3. Evolution of normalised Faraday's law 

The mean electrostatic potential (d>*) is obtained from mean radial momentum balance. Tak- 
ing a flux surface average of the covariant p component of the normalised momentum equation 
[Eq. d45lll. neglecting contributions to the pressure gradient term that are nonlinear in flux surface 
variations of n* and T*, and using the fact that the flux surface average of the radial component 
of Vi must vanish on turbulent timescales to ensure ambipolarity, we obtain 



: 



W* + ^B* 

v A 



X V; 



fro d(p*) p d{<F) l gftp 
(n*) dp dp 2 dp 



Rearranging and integrating with respect to p, we obtain the mean electrostatic potential: 
ip 1 d<n?\ 1 , - 1 rP 



<<*>*> 



eO JO (K) d P 2B eO 1 ' PeO Jo \ { 



Cl>ri • 

W + — B* 

VA 



X V; 



dp. 



In the standard version of CENTORI we define the total electrostatic potential <D* using the adia- 
baticity relation 

/ n * \ 

(60) 



(61) 



which follows from electron force balance along the magnetic field in the limit of vanishing 
electron mass B23I1 . We can write 

dA* _ dA* 

~~dT ~ ~dT' 

where A* is the fluctuating part of A*. It follows from Eqs. (l58l and ( 1591 that 



1 dA* 
va dt 



v* x B* + 



cTeQ (Vp* 



v A Boe\ n\ 

The scalar product of this equation with B* yields 

cT e o 



V(D' 



— B* • — — + B* • E* res 
va dt 



v a Bo e 



B 
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n* 



' £> res ■ 



(62) 



Approximating the time-dependent terms on the left hand side by replacing B* with B* q , we 
obtain 



1 9 (n* a~A cTm 



v A B e I eq \ nl 



We 



Ka ■ I -T 1 - V(& * I + B* • | - VO* |^ - B* q ■ E* res , 



where B* is the normalised fluctuating part of the magnetic field. Neglecting magnetosonic waves 
(i.e. the poloidal component of A*), the covariant representation of A* reduces to 

A* = a;v^. 



In this limit 



and hence 



B* = V x A* = V x (a*V^) = VA* x V£, 



1 d / ~a cT e0 



v A B e P" 



— ^ - VO* + V£ ■ — ^ - VO* x VA* > - B* ■ E* res . 



Using the expression for O* [Eq. (f60b l and the fact that B* q ■ V(/) = for any / since V(/) is 
purely radial and B* q has no radial component, we obtain 



it 



B eq 



B; q -|-^-vo*| = -^-V(n:r;), 



n: 



and hence, using the fact that B* q ■ A* = FA*/(BqR) 



1 F 



dA* 



cT e0 fB, 



^•V(«:r e *) + v^-. 



- VO* x VA* 



va BqR 2 dt va Bq e \ n 
We represent -B* q • E*, es as a diffusion term in this equation by writing 

F 



-B* q -E* res . (63) 



-B; q -E* res - 



RoJD„ 



' d 2 A*. ( Bp ) 2 d 2 A}) 



B 2 do 2 



where the normalised resistive diffusivity is given by 



D„(p,e,o = v* + -\i + 



?W\/^[/yyJ* 2 + W* 2 lV 



X n being a user-defined diffusivity. Turbulent diffusivity terms are present to damp out fluctua- 
tions occurring at the smallest length scales; these tend to be in the poloidal direction, close to 
the magnetic axis. Introducing a parameter Ma = BqR 2 /F we can write 



1 dA* 
v A dt 



cT, 



Ma 



v A B e 
+RoJD„ 



B 



eq 



v«r*) + vf . 



VO* I x VA* 



'd 2 ^ {Bp ) 2 d 2 A*) 
dp 2 + B 2 p dO 2 
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Finally, as in the case of the momentum equation [cf. Eq. d49bl. we add numerical diffusion 
terms to the right hand side of the normalised Faraday's law. 

The above equation is approximated by a finite difference equation which may be written in 
the one-dimensional tridiagonal matrix form 

MjA*™ + Sj A*f 1 + CjA*^ = K Jt (64) 

where subscripts and superscripts have the same meaning as those in the finite difference approx- 
imation to the momentum equation, the coefficients £j and Cj do not depend explicitly on 
A* while %j depends on the latest estimate of this quantity as well as its value at the old timestep 
and also the latest estimate and old value of A* e . As in the case of the velocity components, At at 
the new time is determined by solving the above tridiagonal matrix equation using a predictor- 
corrector scheme. To ensure that only the fluctuating part is actually evolved, the flux surface 
average of A J is evaluated, and subtracted from A* to determine A* at the new time. 

5.3.4. Plasma resistivity 



In terms of the flux surface-averaged density, the Spitzer resistivity is given by [24] 



ill 

^ teer(p) = teH^Z' (65) 

where 

T ce {p) = — — -, (66) 

4 ^(n e )Ae 4 

is the electron collision time, A being the Coulomb logarithm. In a toroidal plasma this is modi- 
fied by neoclassical effects, which, for singly-charged ions, we model using the expression 



(1 _ e l/2)2 + y * 



K v = n^^> (67) 



where _ 

V e = 1/2 ' (68) 

is the dimensionless electron collisionality, q being the safety factor of the flux surface in ques- 
tion, and e — ap 1 ^ 2 /Rq is the local inverse aspect ratio. The resistivity is thus 

7](fi) = ^77Spitzer(p)- (69) 

In the banana regime (v* <K 1) the above expression for K n yields a resistivity which has the 
appropriate limiting behaviour as e — > and e — > 1 12511 : the v* dependence ensures that K n — > 1 
in the limit of high collisionality, as required. 

5.3.5. Evolution of toroidal field parameter F 

The loop voltage Vf is related via the resistive MHD form of Ohm's law to the part of the 
toroidal current associated with FF': 

2V F 

— = -FF'. 

CT] 
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Defining V* F — Vf/2k, it is straightforward to show that the above equation has the following 
solution for F: 



where F vac is the vacuum value of F, i.e. F vac — F at and outside the plasma boundary. The 
plus/minus sign in this expression takes into account the possibility of a reversal in the sign of 
the toroidal magnetic field. 

5.4. Normalised Energy Equations 

We consider here the electron energy equation; the ion equation is treated in a similar manner. 
From Eq. (l39i > we have 

3dT e 3 1 S e 

+ 7v e • vr e + T e V ■ V e = — v ■ q e + — . 

I at I n e n e 

We can write the first term on the right hand side as 

— V ■ q e = -v el] (T e - (T e )) + V • (X e VT e ), 

n e 

where X e is the electron thermal conductivity and v e \\ is the parallel electron thermal relaxation 
rate. The latter may be represented by the expression 

with/ y . a user-defined multiplier and v r / I ( ,(p) = (2(T e )/m e ) 1 ^ 2 the electron thermal velocity. This 
term has the effect of equilibrating the fluctuating component of T e rapidly along the field lines 
at a rate given by v e \\. The portion involving X e is treated as a diffusion term: 

( d 2 T (B ) 2 d 2 T 1 

V ■ (X e VT e ) * R Q J^ e +XRR ) +Xe ^-—tj . 



The Rechester-Rosenbluth diffusivity xrr can be written as [21]: 



B 2 

XRR = fgR Y4 gW n g"' m ' < 71 > 
where < /rr < 1 is a user-defined multiplier. Thus the electron energy equation becomes 

r-jr = ~Ve\\Te+Ve\\(T e )~T e V-\ e --\ e -VT e +—+R JU(Xe +Xrr)^t + ~5T \Xe 



2 dt "«--'w> 2 ' e e n e u " Ann/ dp 2 / B 2 , \ At 50 2 

(72) 

The external source term for this equation is 

S e — P aux.e F"ei T]J » 

where P m x,e is the external heating power per unit volume provided to the electrons (see Sec- 
tion l5.4.2T i. The second term in the expression for S e is the electron-ion equilibration power [Eq. 
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[77} ], through which energy is transferred between the electrons and ions due to the temperature 
difference between them, and the final term is the Ohmic heating power. The external source 
term for the ion energy equation is 

S i — PauxJ + Peii (73) 

where P aux j is the external heating power per unit volume provided to the ions, and P e \ - -P K 
is the electron-ion equilibration power. 

We define the following quantities (with the dimensions of reciprocal length): 

_ V «H r> - % e r> - X RR c * - S 
= > U e = — , Urr = , 0„ = 



v A v a v a v A T e0 n e 

The electron energy equation can then be written in the following normalised form: 
1 r)T* 1 7 1 9 V* 



dp 2 j Bj\ e 86 2 
In a similar fashion we obtain the normalised ion energy equation: 

1 dT* 2 2 2 2 5* 

v.-i = -r; T ' + 3 v i <7 '" > -3 7 '' (V ' v ' +v ' r » ) - v '' V7 '" + 3^ 



3 [dp 2 Bj dO 2 

The Ware pinch term V • T* w , which is only present in the ion equation, is the divergence of the 

flux Hi 

2.44e'/ 2 n'c V F VJ_ 
W v a \ B eq ,pol\ 2nR W 

and the parallel ion thermal relaxation rate \>i\\ is given by 

'"-^{mYik,' <75) 

with fy n a user-defined multiplier. The normalised rate, which again has the dimensions of a 
reciprocal length, is given by vl = v^/va- 

The normalised electron energy equation is approximated by a finite difference equation, with 
the diffusion terms treated exactly by analogy with those in the momentum equation. This can 
be written in tridiagonal matrix form, and solved at each (9, ij) point to advance the normalised 
electron temperature T* at the new time. The normalised ion temperature T* is similarly updated. 

5.4.1. Transport of energy 

The electron collision time is given by Eq. d66i > and the ion collision time by the expression 

3^hT i {T i ? 12 

Tci(p) = a n \>74 4 ' W 
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where Z, is the ion charge state. The power density transferred from electrons to ions (or vice 
versa) due to the temperature difference between them is given by 



Peiip) = — — «T e ) - <r,» = — ^ ^ K) «r;> - <r;» . (77) 

We define dimensionless collisionalities for the two species by the expressions 

yl2{q)R Q <2{q)R 

y e (p) = 175 . v i(P) = -JTo ' (78) 

where v f / v = (2(r,)/m,) 1 ^ 2 is the ion thermal speed. We define flux surface-averaged electron 
and ion cyclotron frequencies and thermal Larmor radii by 

, , e ( B ) , , Zje(B) v t h,e , . v, hJ 

COceKP) = , Uci(P) = . Pe(p) = , PiiP) = • 

m e c m-fi u) ce CD c j 

We also define poloidal Larmor radii by the expressions 

Ppe(p) = Pe — , Ppi(f) = Pi — ■ 

{Hp) \o p ) 
The electron and ion neoclassical thermal diffusivities are taken to be 

K NC ,ee U V pe K NCje ^V 

Xe,N C (p) = -, Xi,N C (p) = , (79) 

Tee ' ci 

where K^a is given by an expression that was proposed by Chang and Hinton 12611 as a finite 
aspect ratio generalisation of a result originally obtained by Hinton and Hazeltine 12711 

_ 0.66+ 1.88e'/ 2 -1.54e 0.66 ((0-74)Vy*) 

Kncj(p) " i+ V^T + o.3iv; + 0311+0.74^/2' 

An identical expression is used for K^c,e, with v* replacing v*. Heat transport in tokamak plasmas 
is typically found to be due mainly to turbulence rather than neoclassical effects, particularly in 
the case of electrons. In MAST ion heat transport can be close to neoclassical in the plasma core 
112811 . where the approximations used to obtain the above expression for K^cj are well-satisfied. 
Closer to the plasma edge in MAST, the ion heat transport is generally dominated by turbulence. 
The thermal diffusivities used in the energy equations have the dimensions of length: 

D e (p, 6,0=^ {^' NC + Xe(^+ q{Rf ^ [fjJ J* 2 + W* 2 ]JJ , (80) 

D i (p,6,0 = -±- {^c+^(l+?W 2 ^[/yyJ* 2 + W* 2 ]JJ, (81) 

where^- e and^-,- are background diffusivities specified by the user. Turbulent diffusivity terms are 
present in Eqs. (l80l and ( ISTT l to damp out fluctuations occurring at suitably small length scales. 
These model phenomenologically the effect of all fluctuations on subgrid scales, in a manner 
similar to that used in large-eddy simulations in meteorology lE9ll . In future work we intend to 
derive suitable closure relations by means of kinetic modelling on scales below those resolvable 
using CENTDRI. 
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5.4.2. Auxiliary Heating Power 

There are three options for the auxiliary electron heating power density profile in CENTORI: 

f (1 -p)P eQe - a Pe^-P^\ 

Pau,,e(p) = j (l-p)PeO B^"***** (82) 

where P e o gives the height of the profile in erg cirT 3 s _1 , a pe is the profile index, and p pea k,e is 
location (~ (r/a) 2 ) at which the power profile peaks. These parameters, along with the choice 
of profile type, are specified by the user. The ion heating profile is treated similarly, with an 
equivalent set of parameters. In principle it is possible to use profiles obtained from radio- 
frequency or neutral beam heating codes (applied to GRASS equilibria), and it is essential to do 
so if precise comparisons with experimental results are required. 

5.5. Normalised mass continuity 

Dividing Eq. d40b by to; TQ va, we obtain the normalised mass continuity equation 

lfti| . . « . . . ,2 17 d 2 n*\ (B p } 2 ( d 2 nl\\ 

(83) 

where the normalised particle source S* (see Section l5.5.U is given by 

S* n =-^- (cm- 1 ). (84) 
m n e v A 

The mass continuity equation is approximated by a finite difference equation which, as in the case 
of the other primary quantities, can be written in a tridiagonal matrix form suitable for advancing 
in time. 

5.5.1. Particle source rate 

The rate at which particles (ions) are supplied externally to the plasma per unit volume is 
5„(p)/m,'. We assume that there are two contributions to this - from an auxiliary (neutral beam) 
power source, if any, and via a density feedback mechanism (see Section I5.5.31 l. The latter 
contribution may be assumed to be highest at the edge, falling to close to zero at the plasma 
centre. Thus, the total normalised particle source rate S* [Eq. d84b l is specified in CENTORI as 

ci*/ \ _ S nip) 1 / Pauxjip) + Paux.eip) c „. . 5fo-l)l /r> C i 

S n ip)= — = — = = — +S nedge C(p)e v >\, (85) 

mjn e v A n e v A \ E beam ' j 

where 5 ne dge is specified in units of cirT 3 s _1 , P fli , A ,,/ e is the external heating power provided to 
the ions/electrons in ergs cirT 3 s _1 , Ebeam is the neutral beam particle energy in ergs, and Cip) is a 
cut-off function used to provide further modulation of the feedback source. Currently we remove 
the feedback source completely outside the p 1 ^ 2 = 0.95 contour, i.e. Cip) — 1 if p 1 ^ 2 < 0.95 and 
Cip) = otherwise. 
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5.5.2. Particle diffusion 

We take the normalised particle diffusivity to be related to the normalised electron thermal 
diffusivity [Eq. <ID]: 

D n ( P , e >0=-^ \xe,NC +Xne[l+ ?<R> 2 ^ [fjj J*' + W* 2 ])J , (86) 

where Xne is a user-defined particle diffusivity. As in the case of Xe and Xi m the thermal dif- 
fusivity expressions [Eqs. (f80T > and (fSTt l. this is used to model transport arising from processes 
occurring on sub-grid scales; typically Xne ~ 10 cm 2 s _1 . 

5.5.3. Density feedback 

There is an option in CENTDRI to use a feedback mechanism to control the volume-averaged 
particle density. This is achieved by modifying the edge particle source rate 5 ne dge at each 
timestep as follows: 

c / (^e target - AW/V) /T OT if «e target > Motal/V 

ned ° e \ otherwise 

where « e target is the requested average density, A^totai is the total number of particles in the plasma 
(i.e. the volume integral of n e ), V is the plasma volume, and r„, is the required timescale for the 
density to reach the target value. If the density is too high the particle source is turned off. 



6. Initial and boundary conditions 

6.1. Initial conditions 

At t = the physical quantities are prescribed as follows. All fluctuating components are 
initialised to zero, except for n e , which is given an arbitrary variation in all three directions. 

v ;,normal(A d > = v u tangential^* d > = 0, v. toroida i(p, 8, = v* v A e'"^, 

n e (p, 6, = n M e- a "P, T e (p, 8, = T M ^ ff " p , T t (p, 8, = T i0 e^, 

The coefficients and profile indices in the above expressions are specified by the user. The ini- 
tial vector potential A and magnetic field B are derived from the initial equilibrium 4r(R, Z), as 
described in Section [2] In the early stages of a simulation it may be necessary to determine 
an equilibrium relatively frequently (typically once every 10 3 time steps) to allow transients to 
settle. This early-stage evolution does not simulate accurately the startup phase of a real plasma. 

6.2. Boundary conditions 

The boundary conditions in the 9 and £ directions are, of course, periodic. In this section we 
discuss the boundary conditions to be applied in the radial direction. 
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6.2.7. Axis boundary conditions 

At each discrete toroidal location the true plasma axis (p = 0, 8, £ = £„) is a coordinate 
singularity, since 8 is undefined (i.e. it can take any value from to 2n). The radial and poloidal 
directions are similarly undefined. There is still a clearly-defined toroidal direction, however. 
With these considerations in mind, the physical components of all vector quantities at the plasma 
axis are dealt with as follows. If V(p, 8, denotes any vector quantity, and p - Ap denotes the 
radial location of the first grid point away from the axis, then the normal and toroidal vector 
components are given by 

y normal(°' ' f») = mean value of ^normal^' ' 

y toroidal(°' °> &) = mean value of y toroidal( A P' ^ &)> 
while the tangential component is set equal to zero. As previously noted, the value of p closest 
to the axis has a small positive value. The scalar quantities n e , T e , Ti, p e and p, are treated in the 
same way as V norma j and Koroidal' wr, il e tne A ux surface-averaged profiles of these quantities 
are assumed to be flat at the magnetic axis. In the case of the density profile, for example, 

<Be>(0) = (n e )(Ap). 

Similar boundary conditions are applied at the axis to (T e ), <r,), (p e ), {p e ) and (p t ). 

6.2.2. Edge boundary conditions 

The edge of the plasma is less problematic in terms of the coordinate system than the axis. 
All four of the following boundary conditions are used for different quantities / in the code: 

• Zero: f(p = 1,8,0 = 0. 

• Flat gradient: df/dp = 0, i.e. f(p = 1, 8, = f(p = 1 - Ap, 8, £). 

• Continuous gradient: df/dp is constant, i.e. d 2 f/dp 2 = 0: 

/(p= 1,8,0 = 2f(p=l-Ap,8,0-f(p=l-2Ap,8,0 (88) 

• Fixed: f(p = 1, 8, is held fixed at some predetermined value. 
These boundary conditions are applied as shown in Table Q] 

7. Evolution of mean and fluctuating components 

7.1. Scalar quantities 

In Section [5] we discussed the equations governing the evolution of physics quantities in 
CENTORI. Each of these quantities can be split into mean (or equilibrium) and fluctuating parts. 
The "mean" of a scalar quantity / in this context simply refers to its flux surface average, as 
defined by Eq. d33t . and the fluctuating component / is the remainder: 

/total = </>+/■ 

In the case of normalised electron density, for example, we have 

n* e (p,8,0 = «)(p) + n* e (p,8,0. 

The normalised quantities T*, T*, p* e and p* are split in a similar fashion. In each case the flux 
surface average is evaluated at each timestep after the total quantity has been updated, and the 
fluctuating component is obtained simply by subtracting the average from the total. 
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quantity 



edge boundary condition 



Af> W ' \ v j,normal' 



zero 



B 
J 

Vi 



(contravariant) continuous gradient 
(contravariant) flat gradient 

(physical) flat gradient (but toroidal component zero) 

(physical) flat gradient 

fixed 

continuous gradient 
continuous gradient 




^.tangential)' ^.toroidal) 



Table 1: Plasma edge boundary conditions applied to evolving quantities in CENTQRI. 



7.2. Vector quantities 

The fluctuating components of vector quantities are obtained in a similar fashion by sub- 
traction of means from totals, but the means themselves are calculated differently. The physical 
components of the mean ion velocity v, e q are given by the flux surface averages of the corre- 
sponding components of the total ion velocity Vi. The mean electron velocity \ e eq , on the other 
hand, is obtained from v, eq and J eq using the flux surface-averaged form of Eq. d55l l. 

The electromagnetic equilibrium vector quantities only need to be re-evaluated when the 
plasma equilibrium is updated (see Section[9]), i.e. when t//(R,Z) is recalculated. Then, the mean 
vector potential A eq is determined using Eqs. ([8j, ( flOb and dT3T >. The mean magnetic field B eq 
is obtained directly from the curl of A eq , and the mean current density J eq is obtained from B eq 
via Ampere's law [Eq. d43tl. However, Eq. ( TT3l > shows that the covariant 6 component of A eq 
depends on F{i[>), which determines the toroidal magnetic field [cf. Eq. ©]. The evolution of F 
is described in Section l5".3.5l 

8. Global energy-related quantities 

The Ohmic heating power density is 



F ohm — "H J — 




The kinetic energy densities in the electrons and ions are given by 




The total thermal energy and magnetic field energy are 




where p is the total pressure. We define the total plasma beta as 



/? = 8tt 



jB 2 dV 3 E B ■ 
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Similarly, the poloidal beta is defined to be 



_ infpdV _ l6nE th 

Pp 



I BjdV 3 j B 2 p dV' 

9. Equilibrium force balance and the Grad-Shafranov equation 

The Grad-Shafranov equation, which can be derived from the steady-state form of the two- 
fluid equations B23I1 . describes the equilibrium state of a current-carrying magnetised plasma 
in which the Lorentz force is balanced by a pressure gradient force. As described below, a 
pseudo-transient approach is used in CENTORI to solve this equation. Similar techniques have 
been employed in computational fluid dynamics [30], but have not, as far we are aware, been 
applied previously to the problem of determining toroidal plasma equilibria. The Grad-Shafranov 
equation can be generalised to include transonic flows and momentum sources IU9I1 . Currently, 
however, only the simplest form of the equation, which is applicable when toroidal flows are 
subsonic and poloidal flows are less than the sound speed multiplied by the ratio of the poloidal 
magnetic field to the total field l3~lll . is used in CENTORI; it can be written in the form 

dR\RdR) dZ 2 F 
where primes denote derivatives with respect to iff. The equation can also be written in the form 

—J ( = -4nR 2 p' - FF', (90) 



c 



where = (c/4n)A*i// is the covariant £ component of the equilibrium current density, J eq . 
Although flow modifications to equilibrium flux surfaces are neglected in the current version of 
CENTORI, the effects of low Mach number flows and flow shear on turbulence and MHD insta- 
bilities are taken into account in the two-fluid equations described in Section|4] Thus, CENTORI 
can be used to model, amongst other things, the stabilising effects of sheared flows on ion tem- 
perature gradient modes 13211 and the destabilising effects of such flows on Kelvin-Helmholtz 
instabilities [33]. It is anticipated that flow effects on plasma equilibria will be taken into ac- 
count in future versions of the code; users of the present version should note that it is strictly 
applicable only to subsonic equilibrium flows. 

9.1. The GRASS free boundary equilibrium solver 

The CENTORI source code includes a free boundary Grad Shafranov equilibrium solver named 



GRASS [34] (GRAd Shafranov Solver), which is used to compute solutions of Eq. ( |89l ), taking 
into account the presence of currents in poloidal field coils. Figure [3] shows the layout of the 
computational domain used in this subroutine. The toroidal field coils are assumed to lie entirely 
outside the computational domain; as described in Section l5.3.51 the toroidal field parameter F 
is determined by the loop voltage and the resistivity. The solver uses two rectangular grids: 

1. The main solution grid, within the domain (R m \ n ,Z m i n ) to (R m ax,Z max ). The plasma and 
the coils are assumed to lie wholly within this grid. The poloidal flux values on the grid 
boundaries t]/- m , (Aout, 'Atop and t^ot are calculated analytically from the given coil currents 
and an approximation to the current distribution in the plasma region. 
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Figure 3: Schematic diagram of the computational domain used in the GRASS equilibrium solver, showing the main solution 
grid and the plasma mask. 



2. The plasma mask, comprising the rectangular region (R pm in,Z pm i n ) to (R pm ax,Z pmax ). The 
mask must not extend outside the main solution grid. The (hot) plasma is assumed to lie 
wholly within the plasma mask, but no coils can be present inside it. 

The coils' current density J c (which needn't be the same in each coil) is assigned to a number of 
grid cells, to approximate the coil locations and cross-section areas. 

9.1.1. GRASS solution procedure 

It is necessary to solve the following equation over the main solution grid: 

A> = — RJ tor , 

c 

where J, or is a function of if/ throughout the region containing plasma, and J lor = J c at the coil 
locations. Thus we can rewrite the equation as 

\n 

A>= —(RJ,H + RJ C ), (91) 
c 



where 



_ , 1 inside plasma mask 
' elsewhere 



and J, is the toroidal component of J eq within the plasma: 



RJ, = -cR 2 p' - ^-FF'. (93) 

4-7T 
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We denote by Ef the toroidal electric field that drives the portion of the toroidal current density 
proportional to —FF'. From the resistive MHD form of Ohm's law we thus have 

4tt , , 4nRE F 
—RJ, = -4nR 2 p' + -. 

C C 7] 

Setting 2nR E F = W, the equivalent loop voltage, we obtain 

cR 2 dp V* 

RJ t = -—JL + -L, (94) 
Ai/f dp i] 

where V* F = Vf/2jt and At/r = tAedge - <Amin (see below). The dependencies of dp/dp and tj on if/ 
are prescribed. 

To determine V* F we divide Eq. (|94l by R and integrate over the poloidal cross-section area, 
identifying this quantity as the total plasma current I p : 

f rJ „ c cRd p^ c v *f JA c r n d p^ t,» r dA 

I B = J,dA = - -dA+ —dA = R— dA + Vl — . 

' J J Mr dp J Rr] MfJ dp J Rrj 

It follows from this that 

F r dA 

J Rij 

It is convenient to introduce a new dependent variable u - i///R 1 ^ 2 satisfying the boundary con- 
ditions 

<Ain(Z) ,„ <Aout(Z) 

"min "max 

"botW - Rl/2 > M top(«) - 1 • 

Ri 

It is also convenient to express if/ as the sum of two terms: which vanishes at Z = Z mm and 
Z = Z max ; and 

7-7 7 -7 



, — — 'nun , — ii[>i.\ — . 
^2 = J l^top + 7 <Abot, 



The quantity i/^i is then equal to i/r - 0- 2 . Equivalently, 
7-7 7 -7 

«2 — ] M top H ] M bot, 

« h 

and u\ = u — u%. Clearly mi vanishes at Z = Z m j n and Z = Z max , making it possible to compute 
this quantity by applying a sine Fourier transform in Z. 
Defining the operator A* by the equation 

1 d 2 u d 2 u 3 



we find that the Grad-Shafranov equation becomes 

d 2 ui d 2 ui 3 _4n(RJ,H 
i - 
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d 2 u x d 2 u x 3 4n(RJ t H ,,, \ 



Since U2 is a prescribed function of Z, the quantity A*M2 only needs to be evaluated once, at the 
beginning of the calculation. Moreover we can set A*«2 = _ 3«2/(4R 2 ), since it is independent of 
R and depends only linearly on Z. 

We approach the problem of solving Eq. (|95l l by considering the parabolic equation 



— = e (A M «i - or) 



(96) 



where a is the right hand side of Eq. d95l l, e is a prescribed pseudo-conductivity (taken to be 
uniform across the poloidal plane) and t is a fictitious, time-like iteration variable (not to be 
confused with the true time, t). The problem of solving Eq. (l95l l is thus equivalent to finding 
"steady-state" solutions of Eq. (t96b . The sine transform of Eq. (l96*l l can be approximated by the 
finite difference equation 



At 



'll-l 



(AR) 2 



(A/?) 2 



7T 2 k 2 



4R 



where the fc^-th sine transform coefficients are denoted by 77], i labels the z'-th grid point in the 
R direction, with grid spacing AR, N labels the pseudo-time variable, and At is the pseudo-time 
step. This equation can be rearranged to give 



ftf ^ 1 1 1 + AT£ 



(AR) 2 



h 2 



3 
4R 2 



Are 



Ate 



(AR) 



2 M li+l 



(AR) 



2 M li-l 



Are Si, 



which is a tridiagonal matrix equation of the form 



Mi C-i + & C ' + Ci fifi'i = («u " Area;) 



where 



^ = - 



Are 



(A/?) 



2 ' 



3i = 1 + Are 



(A^) 2 



/i 2 4/? 2 



ATe 



(AR) 



2 ' 



The tridiagonal matrix equation is straightforward to solve for fi? f+1 ; the inverse sine transform 
of this yields u\ and thereby ifj\. The total flux ^ is recovered by adding ip2> an d tne process is 
repeated until if/ over the grid does not change significantly between pseudo-time steps. 



9.1.2. Plasma current 

In general the plasma current density J, and Aifr change between successive pseudo-time 
steps, and so the the evolution described by Eq. ( f96b is non-linear. It should be noted that the 
dependencies of dp/dp and 77 on iff (or p) are determined externally using CENTORI, rather than 
GRASS. Ideally these functions should vary with if/ in such a way that the residual plasma current 
outside of the chosen edge plasma contour remains negligible. 

9.1.3. Defining the plasma edge 

Once a convergent solution for the equilibrium has been obtained, it is necessary to locate 
the edge of the plasma, which is defined to lie wholly within the plasma mask. By estimating 
\Wiff\ at all grid points within the mask using finite differences, it is straightforward to find all the 
stationary points of iff; these are either X-points (saddle points) or the magnetic axis, which is 
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defined to lie at the global minimum of if/ within the mask. There should be no other stationary 
points of if/ inside the mask, unless some coils have been erroneously located within it. The 
edge of the plasma is then defined after finding a reference if/ m!lx using the following criteria (the 
situation is topologically more complicated in general, but in practice this algorithm suffices): 

• If there are no X-points within the mask, if/ mdx is chosen to be the minimum value of if/ 
along the plasma mask perimeter or the minimum value of if/ at a user-defined set of (R, Z) 
"limiter" points within the mask, whichever is lower. 

• If any X-points are present, if/ mdx is chosen to be either the if/ of the lowest X-point or the 
minimum value of ifi along the inner or outer edges of the mask or the limiter points, if this 
is lower than the if/ of the lowest X-point. 

This ensures that the if/ mdx contour is the largest closed contour within the mask. We then define 
the plasma edge contour i^ e dge to be 

"Aedge = <Aaxi.s + /('/'max ~ ^axisX 

where / = 0.99 when no X-points are present within the mask and / = 0.90 otherwise. This 
has the effect of moving the effective plasma boundary to a contour lying slightly inside the last 
closed flux surface, which is necessary to ensure that the coordinate system described in Section|2] 
does not become strongly distorted near the plasma edge, and enables us to approximate the 
physics equations with central differences without incurring unacceptably large truncation errors. 

Finally, if/ is redefined within the plasma mask so that the plasma edge corresponds to if/ — 0. 
CENTORI is passed only this modified if/(R,Z) within the masked region (thus excluding the 
coils), interpolated onto a grid of the same size (i.e. with the same number of elements) using 
Chebyshev fits in R and Z. The algorithm for determining plasma-based coordinates described 
in Section |2] works extremely well when if/(R,Z) is specified in this way, and almost invariably 
yields a Jacobian J' that closely approximates a flux function as a result. 

9.1.4. Control of the magnetic axis location 

It is sometimes useful to be able to hold the magnetic axis at a specified (R,Z) position. For 
example, up-down asymmetric plasmas are often vertically unstable, and in such cases it may be 
difficult to obtain a convergent solution for the equilibrium using GRASS unless it is possible to 
control the plasma location during the convergence cycle. 

Applying a vertical magnetic field makes it possible to control the radial position of the 
plasma, as follows. Suppose that there is a source of poloidal flux if/sz of the form 

R 2 

4*BZ - l/fBZO T~2 ? 

where if/gzo is a constant and Rq is a measure of the major radius (e.g. the value of R at the centre 
of the computational grid). Then 

„, dif/ B7 opfszo 
ViAbz = -^-e R = 2R— r e R . 
oR Ri 



Since the poloidal magnetic field is V£ x Vif/ it follows that the field due to if/sz is uniform and 
vertical: 

n 2lf/ B Z0 

~ ~~R 2 ~ &Z ' 

"o 
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The Lorentz force J p i as x B z on the plasma arising from a positive plasma current J p i as = 7 p i as 
will be inwards if if/ BZ o > 0. 

Similarly, an externally-provided radial magnetic field affects the plasma's vertical position. 
The radial field due to a poloidal flux of the form 

Z-Zn 

ft BR = ^BRQ ■ 



Z ma x Z m i n 

is 

R _ -&BR0 
R ~ ~M7 -7 

and the Lorentz force on the plasma in this case is downwards if i// B ro > and / p i as > 0. 

We adjust the values of if/ B R0 and tf/ B zo during each GRASS convergence step by comparing the 
latest calculated position of the magnetic axis (/? ax i s , Z ax i s ) with the target location (/?target> Z targe t), 
and applying a correction to the fluxes as follows: 

Z ax i s — Z tar get 
&BR0 -» </fBR0 +f fid ■ 



<Abzo -» fozo + / «Ao 



^axis ^target 
lx target 



where / «: 1 (typically / ~ 0.02) to ensure that the change in the fluxes is not substantial. At 
t = we assume that ifr BR o = ij/ BZ o = 0. The applied corrections should modify the radial and 
vertical fields in such a way that the magnetic axis is pushed towards the target location. 

It is important to note that the magnetic fields associated with these externally-applied poloidal 
flux components are curl-free and hence current-free, i.e. there are no additional current sources 
within the grid implied by their presence. Experimentally, vertical and radial magnetic field per- 
turbations of this type can be introduced by changing the currents in poloidal field coils, although 
such field perturbations are in general non-uniform and therefore the uniform field perturbations 
discussed here are somewhat idealised. It is straightforward to incorporate the additional fluxes 
into GRASS by simply modifying the boundary conditions at the edge of the computational grid, 
at the start of each convergence loop, as follows: 



1/2 



111 v y mm 



mm 



/ cr. / Z-Zp R 

l/r in (Z) + lfr BR0 - — + lABZO- 2 
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10. Outline of code structure 



70.7. Source files 

The CENTORI code is written in standard Fortran 95 throughout, and is contained within 
some 21 source files. The bulk of these contain utility modules and routines to perform specific 
tasks such as I/O, parallel (MPI) communication, error handling, numerical evaluation (Cheby- 
shev/Fourier fitting, and so on) and other customised but standard functionality. To make the 
code as portable as possible, we have avoided the use of external numerical libraries. Those 
areas of the code in which such libraries might improve performance almost all occur in non- 
parallel segments, i.e. are run only by the "global" processor. Since the execution of the code is 
overwhelmingly dominated by periods of parallel execution, serial optimisation through the use 
of specialised libraries is unlikely to confer significant benefits. 

The physics within the code described in this paper is confined to two source files. One of 
these contains all the routines for initialising and evolving the physical fields, sources, sinks and 
so on, and also the routines for calculating plasma coordinates from the tf/(R, Z) grid. The other 
contains the GRASS free boundary equilibrium solver as described in the previous section. 

70.2. Parallelization model 

CENTORI runs in parallel, with each MPI process advancing the physics quantities in an 
allocated three-dimensional subdomain. Aggregation of quantities such as flux surface or volume 
integrals and averages are performed across appropriate sections of the process population via 
specially-written routines. Halo-swapping is necessarily frequent due to the extensive calculation 
of derivatives. For the small grid sizes that are suitable for MAST simulations, there is a tendency 
for the parallelisation to become communication-limited for relatively low process counts. There 
are, however, two different implementations of the key numerical routines available, which are 
optimised for different local grid sizes 13511 . Figure [4] shows speed-up versus process count 
for simulations performed on HECToR at EPCC and HPC-FF at the Jiilich Supercomputing 
Centre when a computational grid of 129 x 65 X 33 and the "eager" implementation of the 
key routines, which perform better for a large number of processes, are used. The HECToR 
results were obtained using the Phase 2b system, based on Cray's XE6 hardware, which provides 
dual socket nodes with 2.1 GHz AMD Opteron 12-core processors and uses Cray's proprietary 
Gemini Interconnect; the PGI compiler was used. The HPC-FF results were obtained using dual 
socket nodes with 2.93 GHz Intel Xeon X5570 quadcore processors and QDR Infiniband switch 
network; for these runs the Intel compiler was used. In this particular case the application shows 
almost linear speed-up for process numbers of up to 128, and continues to display a significant 
speed-up even for 512 processes. For a grid size of 129 x 129 x 129 an almost linear speed-up is 
observed up to at least 512 processes (the largest number of processes used so far). Such a grid 
is larger than is normally used, but might be employed for ITER simulations in the future. 

The computational domain is decomposed across the requested number of MPI processes in 
a standard Cartesian communicator, with periodicity automatically invoked in the two angular 
directions. Physics considerations suggest that the number of grid intervals Ng and N% in the 
radial, poloidal and toroidal directions should be roughly in the ratio 4 : 2 : 1 (a benchmarking 
study has confirmed that this aspect ratio delivers results that are close to optimal i35lo . In order 
to split each dimension into equal-sized portions across processes, the corresponding number of 
grid intervals must be one greater than a power of two, and the number of processes in each 
dimension must be an exact power of two. Thus, the computational grid used to model MAST 
typically has = 129, Ng = 65, N( = 33 over a corresponding grid of 8 x 4 x 4 processes. 
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Figure 4: Execution time speed-up versus process count on HECToR (red with + symbols) and HPC-FF (blue with □ 
symbols), relative to the time taken for a 16-process run, for a computational grid of 129 X 65 X 33. This is the typical grid 
used for MAST simulations. 



Novel techniques are used to optimise the serial execution (on each parallel process) of the 
numerical scheme within CENTDRI. Specially designed derived datatypes employing advanced 
pointer techniques are used, together with lazy evaluation and the option to use strip-mining to 
tailor the code's vectorisation. All the numerical vector operations (scalar and vector products, 
gradient, divergence and curl derivatives), and many pure scalar or combined scalar-vector op- 
erations, are contained within functional black boxes, hiding the implementation details of the 
underlying datatypes and the potential internal conversions between vector representations from 
the physics programmer. This enables a researcher to convert a complicated physics equation to 



a single line of code with ease. Full details are given in [ 35] . File handling is performed in paral- 
lel, with each process writing its own output data files. However an option is under development 
that uses MPI-IO routines to amalgamate the I/O more efficiently 13611 . 

In contrast, the nature of the GRASS two-dimensional equilibrium solution over the entire 
poloidal cross-section of the plasma (and beyond) means that it is best performed in the (R, Z) 
laboratory coordinates on a single process, as is the subsequent construction of the plasma co- 
ordinate system. This impacts only weakly on the code performance, as it is not necessary to 
recalculate the poloidal flux contours tf/(R,Z) on timescales shorter than many microseconds, 
and thus GRASS is called only after many thousands of evolution timesteps; Af ~ 10~ 9 s is the 
typical timestep used. A typical run requires around 600 MB of RAM, and it takes around 12 
hours on 128 MPI processes to simulate 1 ms of plasma evolution in typical tokamak conditions. 

10.3. Additional features within the CENTORI package 

The CENTORI code is best considered as a complete software package, rather than simply a 
collection of source and input files. In addition to its normal role for compilation, the makefile 
includes a number of utility functions that perform tasks such as automatic generation of the code 
documentation, and the creation of a tar file containing the entire source code, its documentation 
and visualisation files, and the input and output files. This has proved to be of great benefit in 
keeping all of the data from a given run together for archival purposes. 
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The source code is self-documenting to a degree, using an included parser program (autodoc) 
to generate html files for each subprogram from specially-formatted comment lines within the 
code. In addition, a full ETgX manual is rigorously maintained to ensure its continued strict 
agreement with the evolving source code (this paper is an abridged version of the full manual). 

A comprehensive visualisation suite has been developed to allow straightforward interpreta- 
tion of the physics output from CENTORI. The program (CentoriScope) is written in the IDL 
language ID7I1 . Work is in progress to bring CENTORI into the EU Integrated Tokamak Modelling 
(ITM JH) framework. 

The code is maintained within a private Subversion repository. Currently, access to the code 
is obtainable only by prior permission from the authors. 



11. Example outputs 

An early version of CENTORI was used to study wave propagation in the vicinity of mag- 
netic X-points; analytical results for the evolution of perturbed wave energy and plasma kinetic 
energy in the ideal MHD limit were recovered numerically [39]. In this section we present two 
illustrative examples of calculations that can be performed using the full version of the code. 

77.7. Tearing mode in large aspect ratio tokamak plasma 

We demonstrate the capability of using CENTORI to model MHD instabilities by considering 
the example of a tearing mode in a very large aspect ratio (minor radius a = 0.36 m, major radius 
R = 16.8 m), circular cross-section tokamak plasma with a toroidal magnetic field of 9.7 T. 
For this purpose we prescribe an equilibrium with uniform density (10 24 irT 3 ) and temperature 
(T e = Tj = 46 eV). The resistivity, which we take to be given by the Spitzer expression [Eqs. d65l l 
and d66ll1. is then equal to 3.67 x 1CT 16 s, and the Lundquist number S = Anav^Kc 1 !]) ^ 2 x 10 4 . 
The number of radial grid points (256) was chosen to be sufficiently large that the resistive layer 
width d ~ aS~ 2 l 5 - 0.7 cm was well-resolved. The quantity FF', which is proportional to 
the toroidal current density in this large aspect ratio limit, was prescribed to have the following 
profile: 

(1 + sinh 2 [2.09p]} 

where FF'(0) is a constant, chosen to ensure that the corresponding ^-profile remained above 
unity across the plasma, with q — 2 at a normalised minor radius of about 0.65. This configura- 
tion is expected to be unstable to the growth of a tearing mode with poloidal and toroidal mode 
numbers m — 2, n — 1 [40]. For the purpose of this calculation only the generalised Ohm's law 
and the ion momentum equation were used [Eqs. (l36l l and d37V l. Ohm's law was reduced to the 
resistive MHD form, and viscosity was neglected in the momentum equation (except for the nu- 
merical viscosity described in Section l5.1.1l which is required to suppress numerical oscillations, 
but is set at a level which is sufficiently low for the system to be effectively inviscid). The electro- 
static potential in this case was calculated not using Eq. ( f60b but by evolving the perpendicular 
ion velocity, identifying this as an E x B drift, and integrating to obtain O. Nonlinear terms were 
omitted from both Ohm's law and the momentum equation. Modes with m/n equal to values of 
q inside the plasma other than 2, such as the 3/2 mode, can also be unstable in the presence of a 
current density gradient; in general these modes cannot be exluded from simulations performed 
using a non-spectral code such as CENTORI. For this particular simulation, a Fourier filter was 
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therefore applied at the end of each time step to exclude all harmonics other than the dominant 
2/1 mode. 

As expected, the configuration described above was found using CENTORI to be unstable to 
the growth of a 2/1 tearing mode. Figure [5] shows snapshots of Af and O during the tearing 
mode growth. The profiles of these quantities resemble those obtained using CUTIE in a similar 
(although not identical) parameter regime; cf. Fig. 2 in Ref. [41], which was obtained with S = 
2 x 10 4 (defined in terms of the local resistivity at the magnetic axis) and relatively low viscosity 
(it should be noted that a precise comparison between tearing mode calculations performed using 
CENTORI and the CUTIE results presented in Ref. J4jl] is not in fact possible, since in the latter 
case a low aspect ratio (R/a = 2.5) was assumed for the purpose of calculating the ^-profile but 
the flux surfaces, as in all CUTIE simulations, were taken to be concentric circles; in a toroidal 
code such as CENTORI the circular flux surface limit can only be approached by taking the aspect 
ratio to be very large). The growth rate of the mode shown in Fig.[5]is approximately 2.5x 10~ 3 /-za 
where ta — a/vA is the Alfven time; this is comparable to the rate deduced analytically in Ref. 

i.e. jta ~ S 3/5 . It is somewhat lower than the rate found using CUTIE in the low viscosity 
limit with S = 2 x 10 4 at the magnetic axis (yr A - 1.8 x 10~ 2 ) |4l|], but in this calculation the 
local resistivity at the q — 2 surface was higher than the value at the magnetic axis, implying a 
higher 2/1 tearing mode growth rate. 

fluotootino. A, covzela (Wo) t = 3.250E+00 ms IttetO = O.OOOE+00 zeta = O.OOOE+00 total senior potential (Vj t = 3.25DE+00 no tteto = O.OOOE+00 zeta = O.000E+00 




Figure 5: Radial profiles of Ar (left) and <t> (right) in CENTORI simulation of tearing mode in large aspect ratio tokamak. 



11.2. Turbulence simulation in conventional aspect ratio tokamak plasma 

We present here the results of a CENTORI run simulating 1 ms of a conventional aspect ratio 
tokamak plasma with minor radius 0.55 m, major radius 1 .67 m, elongation 1 .7 and triangularity 
0.18; the equilibrium flux surface contours, computed using GRASS, are shown in Fig. [6] The 
chosen plasma current was 1 MA, the toroidal magnetic field 2.5 T and the plasma volume 15m 3 . 
The initial flux surface-averaged density, temperature and current density (primary quantities) 
were held approximately constant during the simulation by using adaptive sources of the form 
S = -a{(f) - (/o)), where / is the primary variable profile, /b is its initial profile, and a is an 
inverse reaction time response, set equal to the reciprocal of At, the CENTORI time step (0.5 ns, 
in this case). Both the initial ion velocity and the external momentum source S v were set equal 
to zero. The profiles of electron density, electron and ion temperatures, are shown in Fig. [7] 
together with the ^-profile. 

The evolution of the sources follows that of the fluctuations; after an initial transient, lasting 
around lOOps, they reached a quasi-steady level. The boundary conditions were those listed in 
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Figure 6: Plot of if/ contours for plasma equilibrium used in turbulence simulation described in Section 1 1 . The inner rectangle 
indicates the plasma mask used to construct this equilibrium; note that the boundary of the mask lies outside the region of 
confined plasma, bounded by a thick black curve. 



Table 1. The spatial grid comprised 129 radial points, 65 poloidal points and 33 toroidal points. 
The run was executed on 64 processors of the HPC-FF machine at the Jiilich Supercomputing 
Centre, the total wall-clock time being approximately 18 hours. Figure [8] shows the evolution of 
fluctuations in toroidal current density and electron density at p — 0.46, 9 = 0, f = 0. It is evident 
from a comparison of the relative amplitudes of the temporal variations in these two quantities 
that the fluctuations are electromagnetic in character. 

Figure [9] shows a poloidal cross-section of the toroidal current density fluctuations at £ = 0. 
The X symbol in this figure marks the location chosen for the sample of local fluctuations shown 
in Fig. [8] The temporal evolution of the electron thermal conductivity is shown in Fig. |T0j In 
the plasma turbulence literature this quantity is often normalised to the gyro-Bohm diffusivity 
Xgb = P^cJLj where c s = (T e /mi) 1 ^ 2 , p s = c s /u> C j and Lj = T e l{dT e ldr) is the temperature 
scale length 14211 . In the case of the local plasma parameters corresponding to the results shown 
in Fig. [10] p 2 s c s /Lr — 12 m 2 s _1 ; normalised to this value, the time-averaged thermal conductivity 
plotted in Fig. [10] is around 2, which is comparable to normalised Xe values in gyro-kinetic 
simulations reported by Peeters and co-workers B42I1 . 

The results presented above can be used to estimate the magnitudes of the potential and in- 
ductive contributions to the turbulent electric field; as noted in Section 15.1.1 1 only the potential 
electric field term is retained in the ion momentum equation in CENT0RI. From Fig. [8] we note 
that the electron density fluctuations have a relative amplitude h e /n e of the order of 10~ 2 . Elec- 
tron force balance implies that the associated electrostatic potential fluctuations <J) are of order 
10~ 2 r f /e ~ 20 V, since the electron temperature at this point in the plasma is about 2keV (cf. 
Fig- 13 ■ Figure [9] indicates that the fluctuations have a characteristic scale length perpendicu- 
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Figure 7: Profiles used in turbulence simulation described in Section 11. The thick solid curve shows the electron density in 
units of 10"itT 3 , the thin solid curve is the (/-profile, and the dashed and dashed-dotted curves show respectively the electron 
temperature and ion temperature in keV. 



lar to the magnetic field L± of order l(T 2 m, suggesting potential electric field fluctuations of 
<t>/Lj_ ~ 2kVirT 1 . In contrast, the frequency (to ~ 200krad s _1 ) and amplitude (~ 0.03 MArrT 2 ) 
of the current fluctuations / shown in the upper frame of Fig. [8] imply inductive electric fields of 
order uj^Jl}^ ~ 1 VirT 1 Qiq being the permeability of free space). Thus, for the parameters of 
this simulation (which are fairly representative of hot tokamak plasmas), the potential compo- 
nent of the fluctuating electric field is around three orders of magnitude larger than the inductive 
component, and our neglect of the latter in the ion momentum equation is therefore fully justified. 

12. Conclusion 

We have presented a comprehensive description of a novel two-fluid electromagnetic plasma 
turbulence code, CENTDRI, together with sample output from a CENTDRI simulation of a large 
aspect ratio tokamak plasma. The code is used to compute self-consistently the time evolution of 
plasma fluid quantities and fields in a toroidal configuration of arbitrary aspect ratio and plasma 
beta. The code is parallelised, and the equations are represented in fully finite difference form, 
ensuring good scalability. The equations are solved in a plasma coordinate system that is defined 
such that the Jacobian of the transformation from laboratory coordinates is a function only of 
the equilibrium poloidal flux, thereby accelerating vector operations and the evaluation of flux 
surface averages. GRASS, a subroutine of CENTDRI, is used to determine the plasma equilibrium 
(and hence the plasma coordinates in which the fluid and Maxwell equations are solved) by 
computing the steady-state solutions of a diffusion equation with a pseudo-time derivative. The 
physics model implemented in CENTDRI is based solidly on that used in the highly-successful 
CUTIE code, and we are confident that it will prove to be a powerful tool for the study of heat, 
particle and momentum transport in tokamak plasmas. In a forthcoming paper we will report the 
results of the first global simulations, performed using CENT0RI, of electromagnetic, nonlinearly- 
saturated turbulence and transport in a spherical tokamak plasma (MAST). 
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Figure 8: Fluctuations in toroidal current density and electron density in outer midplane at p = 0.46. 
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